! Copyright (C) 2009-2019 University of Warwick
!
! This program is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this program.  If not, see <http://www.gnu.org/licenses/>.

MODULE fields

   USE boundary
   !debug cgy
   USE shared_data
   !GuanyueChen vcm=====
   USE evaluator

   IMPLICIT NONE

   INTEGER :: field_order
   REAL(num) :: hdt, fac
   REAL(num) :: hdtx, hdty
   REAL(num) :: cnx, cny
   REAL(num) :: alphax = 1.0_num, alphay = 1.0_num
   REAL(num) :: betaxy = 0.0_num, betayx = 0.0_num
   REAL(num) :: deltax = 0.0_num, deltay = 0.0_num

   !GuanyueChen vcm=====
   REAL(num) :: exmu_max, eymu_max,ezmu_max
   REAL(num) :: emu_max
   REAL(num) :: emu_max_global
   INTEGER :: step_limit_global
   !REAL(num_hp) :: exmu_max, eymu_max,ezmu_max
   !REAL(num_hp) :: emu_max
   !REAL(num_hp) :: emu_max_global
CONTAINS

!===============main process==============
   SUBROUTINE update_eb_fields_half

      hdt  = 0.5_num * dt
      hdtx = hdt / dx
      hdty = hdt / dy

      cnx = hdtx * c**2
      cny = hdty * c**2

      fac = hdt / epsilon0

      ! Update E field to t+dt/2
      CALL update_e_field

      ! Now have E(t+dt/2), do boundary conditions on E
      CALL efield_bcs

      ! Update B field to t+dt/2 using E(t+dt/2)
      CALL update_b_field

      ! Now have B field at t+dt/2. Do boundary conditions on B
      CALL bfield_bcs(.TRUE.)

      ! Now have E&B fields at t = t+dt/2
      ! Move to particle pusher

   END SUBROUTINE update_eb_fields_half



   SUBROUTINE update_eb_fields_final

      hdt  = 0.5_num * dt
      hdtx = hdt / dx
      hdty = hdt / dy

      cnx = hdtx * c**2
      cny = hdty * c**2

      fac = hdt / epsilon0

      CALL update_b_field

      CALL bfield_final_bcs

      CALL update_e_field

      CALL efield_bcs

   END SUBROUTINE update_eb_fields_final


!===============main vcm process==============
!========GuanyueChen vcm==================
!Added by Jing.W in 2021.1
   SUBROUTINE update_eb_fields_half_qed
      INTEGER :: ix,iy,counter,step_limit
      INTEGER :: handle_block
      hdt  = 0.5_num * dt
      hdtx = hdt / dx
      hdty = hdt / dy

      cnx = hdtx * c**2
      cny = hdty * c**2

      fac = hdt / epsilon0

      !write(*,*) 'ey', ey(-2) ,ey(-1) ,ey(0) ,ey(1), ey(2)
      !write(*,*) 'bz', bz(-2) ,bz(-1) ,bz(0) ,bz(1), bz(2)

      !calculate polarization and magnetization P(n),M(n),qPx0,qPy0,qPz0,qMx0,qMy0,qMz0
      CALL interpolation_for_ex
      CALL calculate_PM_for_ex
      !     write(*,*) 'ey',ey(-1) ,ey(0) ,ey(1), ey(nx-1),ey(nx),ey(nx+1)
      !     write(*,*) 'bz',bz(-1) ,bz(0) ,bz(1), bz(nx-1),bz(nx),bz(nx+1)
      ! write(*,*) 'qPy0', qPy0(0) ,qPy0(1), qPy0(nx), qPy0(nx-1)
      ! write(*,*) 'qMz0', qMz0(0) ,qMz0(1), qMz0(nx), qMz0(nx-1)
      ! Update E field of zero order to t+dt/2 e0, b0
      CALL update_e_field_zero
      ! !GuanyueChen=====
      ! !debug
      ! if (rank==0) then
      !   write(*,*)"运行到第",step,"步，  "
      !   write(*,*)"在rank：",rank
      !   write(*,*)"ey:",bz(:,:)
      !   CALL MPI_SEND(testflag, 1, MPI_INTEGER, rank+1, tag, comm, errcode)
      ! else
      !   CALL MPI_RECV(debug_temp, 1, MPI_INTEGER, rank-1, tag, comm, status, errcode)
      !   write(*,*)"运行到第",step,"步，  "
      !   write(*,*)"在rank：",rank
      !   write(*,*)"ey:",bz(:,:)
      !   if (rank<(nproc-1)) then
      !     CALL MPI_SEND(testflag, 1, MPI_INTEGER, rank+1, tag, comm, errcode)
      !   endif
      ! endif
      ! if (step==testflag)then
      !   pause
      ! endif

      CALL efield0_bcs

      ! Update B field of zero order to t+dt/2 e0, b0
      CALL update_b_field_zero

      CALL bfield0_bcs(.FALSE.)

      emu_max = 0.0_num
      !emu_max = 0.0_num_hp
      !emu_max_global = 10000.0_num_hp
      emu_max_global = 10000.0_num
      step_limit = 0
      DO while( step_limit  < 100 )
         !P(n+1/2), M(n+1/2)
         step_limit = step_limit+1
         CALL interpolation_for_ex0
         CALL calculate_PM_for_ex0

         !write(*,*) 'qPy', step_limit
         ! if(step_limit == 1) then
         ! write(*,*) 'ey0',ey0(-1) ,ey0(0) ,ey0(1), ey0(nx-1),ey0(nx),ey0(nx+1)
         ! write(*,*) 'bz0',bz0(-1) ,bz0(0) ,bz0(1), bz0(nx-1),bz0(nx),bz0(nx+1)
         ! write(*,*) 'qPy', qPy(0) ,qPy(1), qPy(nx), qPy(nx-1)
         ! write(*,*) 'qMz', qMz(0) ,qMz(1), qMz(nx), qMz(nx-1)
         ! endif
         !qed-corrected maxwell slover,E(n+1/2),B(n+1/2)--ex1,ey1,ez1,bx1,by1,bz1
         CALL update_e_field_vcm
         CALL efield1_bcs
         CALL update_b_field_vcm
         CALL bfield1_bcs(.FALSE.)

         exmu_max = MAXVAL(ABS(ex1-ex0))
         eymu_max = MAXVAL(ABS(ey1-ey0))
         ezmu_max = MAXVAL(ABS(ez1-ez0))
         emu_max = MAX(exmu_max,eymu_max,ezmu_max)
         CALL MPI_Allreduce(emu_max, emu_max_global, 1, mpireal, MPI_MAX, comm, errcode)

         IF( emu_max_global  < 0.001_num ) EXIT
         DO ix = 1-ng,nx+ng
            do iy = 1-ng, ny+ng
               ex0(ix,iy) = ex1(ix,iy)
               ey0(ix,iy) = ey1(ix,iy)
               ez0(ix,iy) = ez1(ix,iy)
               bx0(ix,iy) = bx1(ix,iy)
               by0(ix,iy) = by1(ix,iy)
               bz0(ix,iy) = bz1(ix,iy)

            end do

         ENDDO
      ENDDO

      CALL MPI_Allreduce(step_limit, step_limit_global, 1, MPI_INTEGER , MPI_MAX, comm, errcode)
      if (x_max_boundary .and. y_max_boundary ) then
         write(*,*)'in update_eb_fields_half_qed',step_limit_global
         write(*,*)'final max global error:',emu_max_global
      endif

      DO ix = 1,nx
         do iy = 1, ny
            ex(ix,iy) = ex1(ix,iy)
            ey(ix,iy) = ey1(ix,iy)
            ez(ix,iy) = ez1(ix,iy)
         end do
      ENDDO
      !write(*,*) 'ey', ey(-2) ,ey(-1) ,ey(0) ,ey(1), ey(2)
      !write(*,*) 'bz', bz(-2) ,bz(-1) ,bz(0) ,bz(1), bz(2)
      ! Now have E(t+dt/2), do boundary conditions on E
      CALL efield_bcs
      ! write(*,*) 'ey',ey(-1) ,ey(0) ,ey(1), ey(nx-1),ey(nx),ey(nx+1)

      ! Update B field to t+dt/2 using E(t+dt/2)
      CALL update_b_field

      ! Now have B field at t+dt/2. Do boundary conditions on B
      CALL bfield_bcs(.TRUE.)

      !write(*,*) 'ey', ey(-2) ,ey(-1) ,ey(0) ,ey(1), ey(2)
      !write(*,*) 'bz', bz(-2) ,bz(-1) ,bz(0) ,bz(1), bz(2)

      !write(*,*) bz(0)
      ! Now have E&B fields at t = t+dt/2
      ! Move to particle pusher
   END SUBROUTINE update_eb_fields_half_qed


   SUBROUTINE update_eb_fields_final_qed
      INTEGER :: ix,iy,counter,step_limit
      INTEGER, DIMENSION(1) :: loc
      INTEGER :: handle_block
      hdt  = 0.5_num * dt
      hdtx = hdt / dx
      hdty = hdt / dy

      cnx = hdtx * c**2
      cny = hdty * c**2

      fac = hdt / epsilon0

      DO ix = 1,nx
         do iy = 1, ny
            ex(ix,iy) = ex1(ix,iy)
            ey(ix,iy) = ey1(ix,iy)
            ez(ix,iy) = ez1(ix,iy)
         end do
      ENDDO

      !calculate polarization and magnetization P(n+1/2),M(n+1/2),qP0,qM0
      CALL interpolation_for_ex
      CALL calculate_PM_for_ex
      CALL update_b_field

      CALL bfield_final_bcs

      !update_b_field_zero
      Do ix = 1-ng,nx+ng
         DO iy = 1-ng,ny+ng
            bx0(ix,iy) = bx(ix,iy)
            by0(ix,iy) = by(ix,iy)
            bz0(ix,iy) = bz(ix,iy)
         enddo
         !write(*,*) ix
      ENDDO
      !Updata E field of zero order to t+dt, e0, b0
      CALL update_e_field_zero
      !calculate boundary for e0!!!!!
      CALL efield0_bcs
      !emu_max = 0.0_num_hp
      emu_max = 0.0_num
      !emu_max_global = 10000.0_num_hp
      emu_max_global = 10000.0_num
      step_limit = 0
      DO while( step_limit  < 100 )
         step_limit = step_limit+1
         !using B(n+1) and E_0(n+1) --> P(n+1), M(n+1)
         CALL interpolation_for_ex0
         CALL calculate_PM_for_ex0
         !qed-corrected maxwell slover
         CALL update_e_field_vcm
         !calculate boundary for e1 !!!!!
         CALL efield1_bcs

         exmu_max = MAXVAL(ABS(ex1-ex0))
         eymu_max = MAXVAL(ABS(ey1-ey0))
         ezmu_max = MAXVAL(ABS(ez1-ez0))

         if (step_limit > 48 ) THEN
            loc = MAXLOC(ABS(ey1-ey0),1)
            !write(*,*)   step_limit
            !write(*,*)  'delta_ey',  delta_ey(100)
            !write(*,*)  'delta_ey',  delta_ey(loc(1)-3)
            !write(*,*)  'ey',  ey(100)
            !write(*,*)  'ey+delta_ey',  ey(loc(1)-3)+delta_ey(loc(1)-3)
            !write(*,*)  'loc',  loc(1)
            !write(*,*)  'ey1',  ey1(loc(1)-3)
            !write(*,*)  'ey0',  ey0(loc(1)-3)
            !write(*,*)  'eymu', ey1(loc(1)-3)-ey0(loc(1)-3)
            !write(*,*) 'emu_max_global', eymu_max
            !write(*,*)  MAXLOC(ABS(ey1-ey0))
         ENDIF
         emu_max = MAX(exmu_max,eymu_max,ezmu_max)
         CALL MPI_Allreduce(emu_max, emu_max_global, 1, mpireal, MPI_MAX, comm, errcode)

         !write(*,*) 'emu_max_global', emu_max_global
         IF( emu_max_global  < 0.001_num ) EXIT



         !emu_max0 = emu_max
         DO ix = 1-ng,nx+ng
            DO iy = 1-ng,ny+ng
               ex0(ix,iy) = ex1(ix,iy)
               ey0(ix,iy) = ey1(ix,iy)
               ez0(ix,iy) = ez1(ix,iy)
            ENDDO
         ENDDO
         !write(*,*)  ex0(200)
         !write  ex0(200)
      ENDDO
      CALL MPI_Allreduce(step_limit, step_limit_global, 1, MPI_INTEGER , MPI_MAX, comm, errcode)
      if (x_max_boundary .and. y_max_boundary ) then
         write(*,*)'in update_eb_fields_final_qed',step_limit_global
         write(*,*)'final max global error:',emu_max_global
      endif

      DO ix = 1-ng,nx+ng
         DO iy = 1-ng,ny+ng
            ex(ix,iy) = ex1(ix,iy)
            ey(ix,iy) = ey1(ix,iy)
            ez(ix,iy) = ez1(ix,iy)
         ENDDO
      ENDDO

      CALL efield_bcs

   END SUBROUTINE update_eb_fields_final_qed


!=============vcm field advance=========
   SUBROUTINE interpolation_for_ex
      TYPE(parameter_pack) :: parameters
      INTEGER :: ix,iy,err_livky=0
      DO ix = 1-ng, nx+ng
         DO iy= 1-ng,ny+ng
            if (static_handle==2  )then
               parameters%pack_iy = iy
               parameters%pack_ix = ix
               static_ey(ix,iy)=evaluate_with_parameters(static_ey_stack, parameters, err_livky)
               static_bz(ix,iy)=evaluate_with_parameters(static_bz_stack, parameters, err_livky)
            endif
            qex(2*ix,2*iy) = ex(ix,iy)+static_ex(ix,iy)
            qey(2*ix,2*iy) = ey(ix,iy)+static_ey(ix,iy)
            qez(2*ix,2*iy) = ez(ix,iy)+static_ez(ix,iy)
            qbx(2*ix,2*iy) = bx(ix,iy)+static_bx(ix,iy)
            qby(2*ix,2*iy) = by(ix,iy)+static_by(ix,iy)
            qbz(2*ix,2*iy) = bz(ix,iy)+static_bz(ix,iy)
         ENDDO
      ENDDO

      DO ix = 1-ng, nx+ng-1
         DO iy= 1-ng,ny+ng-1
            qex(2*ix+1,2*iy) = 0.5_num * (ex(ix,iy)+ex(ix+1,iy))+0.5_num * (static_ex(ix,iy)+ static_ex(ix+1,iy))
            qey(2*ix+1,2*iy) = 0.5_num * (ey(ix,iy)+ey(ix+1,iy))+0.5_num * (static_ey(ix,iy)+ static_ey(ix+1,iy))
            qez(2*ix+1,2*iy) = 0.5_num * (ez(ix,iy)+ez(ix+1,iy))+0.5_num * (static_ez(ix,iy)+ static_ez(ix+1,iy))
            qbx(2*ix+1,2*iy) = 0.5_num * (bx(ix,iy)+bx(ix+1,iy))+0.5_num * (static_bx(ix,iy)+ static_bx(ix+1,iy))
            qby(2*ix+1,2*iy) = 0.5_num * (by(ix,iy)+by(ix+1,iy))+0.5_num * (static_by(ix,iy)+ static_by(ix+1,iy))
            qbz(2*ix+1,2*iy) = 0.5_num * (bz(ix,iy)+bz(ix+1,iy))+0.5_num * (static_bz(ix,iy)+ static_bz(ix+1,iy))

            qex(2*ix,2*iy+1) = 0.5_num * (ex(ix,iy)+ex(ix,iy+1))+0.5_num * (static_ex(ix,iy)+ static_ex(ix,iy+1))
            qey(2*ix,2*iy+1) = 0.5_num * (ey(ix,iy)+ey(ix,iy+1))+0.5_num * (static_ey(ix,iy)+ static_ey(ix,iy+1))
            qez(2*ix,2*iy+1) = 0.5_num * (ez(ix,iy)+ez(ix,iy+1))+0.5_num * (static_ez(ix,iy)+ static_ez(ix,iy+1))
            qbx(2*ix,2*iy+1) = 0.5_num * (bx(ix,iy)+bx(ix,iy+1))+0.5_num * (static_bx(ix,iy)+ static_bx(ix,iy+1))
            qby(2*ix,2*iy+1) = 0.5_num * (by(ix,iy)+by(ix,iy+1))+0.5_num * (static_by(ix,iy)+ static_by(ix,iy+1))
            qbz(2*ix,2*iy+1) = 0.5_num * (bz(ix,iy)+bz(ix,iy+1))+0.5_num * (static_bz(ix,iy)+ static_bz(ix,iy+1))

            !下面一块代码改为和师姐写法一致，本质相同
            ! qex(2*ix+1,2*iy+1)=0.25_num*(ex(ix,iy)+ex(ix,iy+1)+ex(ix+1,iy)+ex(ix+1,iy+1))
            ! qey(2*ix+1,2*iy+1)=0.25_num*(ey(ix,iy)+ey(ix,iy+1)+ey(ix+1,iy)+ey(ix+1,iy+1))
            ! qez(2*ix+1,2*iy+1)=0.25_num*(ez(ix,iy)+ez(ix,iy+1)+ez(ix+1,iy)+ez(ix+1,iy+1))
            ! qbx(2*ix+1,2*iy+1)=0.25_num*(bx(ix,iy)+bx(ix,iy+1)+bx(ix+1,iy)+bx(ix+1,iy+1))
            ! qby(2*ix+1,2*iy+1)=0.25_num*(by(ix,iy)+by(ix,iy+1)+by(ix+1,iy)+by(ix+1,iy+1))
            ! qbz(2*ix+1,2*iy+1)=0.25_num*(bz(ix,iy)+bz(ix,iy+1)+bz(ix+1,iy)+bz(ix+1,iy+1))
            

         ENDDO
      ENDDO

      Do ix = 1-ng, nx+ng-1
         Do iy = 1-ng, ny+ng-1
             qex(2*ix+1,2*iy+1) =0.25_num * (qex(2*ix,2*iy+1)+qex(2*ix+2,2*iy+1)+qex(2*ix+1,2*iy)+qex(2*ix+1,2*iy+2))
             qey(2*ix+1,2*iy+1) =0.25_num * (qey(2*ix,2*iy+1)+qey(2*ix+2,2*iy+1)+qey(2*ix+1,2*iy)+qey(2*ix+1,2*iy+2))
             qez(2*ix+1,2*iy+1) =0.25_num * (qez(2*ix,2*iy+1)+qez(2*ix+2,2*iy+1)+qez(2*ix+1,2*iy)+qez(2*ix+1,2*iy+2))
             qbx(2*ix+1,2*iy+1) =0.25_num * (qbx(2*ix,2*iy+1)+qbx(2*ix+2,2*iy+1)+qbx(2*ix+1,2*iy)+qbx(2*ix+1,2*iy+2))
             qby(2*ix+1,2*iy+1) =0.25_num * (qby(2*ix,2*iy+1)+qby(2*ix+2,2*iy+1)+qby(2*ix+1,2*iy)+qby(2*ix+1,2*iy+2))
             qbz(2*ix+1,2*iy+1) =0.25_num * (qbz(2*ix,2*iy+1)+qbz(2*ix+2,2*iy+1)+qbz(2*ix+1,2*iy)+qbz(2*ix+1,2*iy+2))
         ENDDO
     ENDDO
   END SUBROUTINE interpolation_for_ex


   SUBROUTINE calculate_PM_for_ex
      INTEGER :: ix,iy
      REAL(num) :: xi
      REAL(num) :: Ex, Ey, Ez, Bx, By, Bz
      !xi = 2*alpha**2*epsilon0**2*h_bar**3/45/m0**4/c**5*1e8/2.285
      !xi = 5.71074693e-45_num
      xi = ksi
      !write(*,*) xi
      DO ix = 1-ng+1, nx+ng-1
         DO iy = 1-ng+1, ny+ng-1
            !for Py, Mx
            !沿x方向做平均
            Ex = 0.5*(qex(2*ix,2*iy+1)+qex(2*ix-2,2*iy+1))
            Ey = qey(2*ix,2*iy)
            Ez = qez(2*ix,2*iy+1)
            Bx = qbx(2*ix,2*iy)
            By = 0.5*(qby(2*ix,2*iy+1)+qby(2*ix-2,2*iy+1))
            Bz = qbz(2*ix-1,2*iy)
            ! write(*,*) Ex,Ey,Ez,Bx,By,Bz
            qPy(ix,iy) = 2*xi*(2*(Ex**2+Ey**2+Ez**2-c**2*(Bx**2+By**2+Bz**2))*Ey + 7*c**2*(Ex*Bx+Ey*By+Ez*Bz)*By)
            qMx(ix,iy) = (-1)*2*xi*c**2*(2*(Ex**2+Ey**2+Ez**2-c**2*(Bx**2+By**2+Bz**2))*Bx-7*(Ex*Bx+Ey*By+Ez*Bz)*Ex)
            !write(*,*) 2*xi*(2*(Ex**2+Ey**2+Ez**2-c**2*(Bx**2+By**2+Bz**2))*Ey)
            !for Pz
            Ex = qex(2*ix-1,2*iy)
            Ey = qey(2*ix,2*iy-1)
            Ez = qez(2*ix,2*iy)
            Bx = qbx(2*ix,2*iy-1)
            By = qby(2*ix-1,2*iy)
            Bz = 0.5*(qbz(2*ix,2*iy-1)+qbz(2*ix-2,2*iy-1))
            ! write(*,*) Ex,Ey,Ez,Bx,By,Bz
            qPz(ix,iy) = 2*xi*(2*(Ex**2+Ey**2+Ez**2-c**2*(Bx**2+By**2+Bz**2))*Ez + 7*c**2*(Ex*Bx+Ey*By+Ez*Bz)*Bz)
            !for Px, My
            Ex = qex(2*ix,2*iy)
            Ey = 0.5*(qey(2*ix,2*iy-1)+qey(2*ix+2,2*iy-1))
            Ez = qez(2*ix+1,2*iy)
            Bx = 0.5*(qbx(2*ix,2*iy-1)+qbx(2*ix+2,2*iy-1))
            By = qby(2*ix,2*iy)
            Bz = qbz(2*ix,2*iy-1)
            qPx(ix,iy) = 2*xi*(2*(Ex**2+Ey**2+Ez**2-c**2*(Bx**2+By**2+Bz**2))*Ex + 7*c**2*(Ex*Bx+Ey*By+Ez*Bz)*Bx)
            qMy(ix,iy) = (-1)*2*xi*c**2*(2*(Ex**2+Ey**2+Ez**2-c**2*(Bx**2+By**2+Bz**2))*By-7*(Ex*Bx+Ey*By+Ez*Bz)*Ey)
            ! write(*,*) Ex,Ey,Ez,Bx,By,Bz
            !write(*,*) qPx0(ix),qPy0(ix),qPz0(ix)
            !write(*,*) qMx0(ix),qMy0(ix),qMz0(ix)
            !for Mz
            Ex = qex(2*ix,2*iy+1)
            Ey = qey(2*ix+1,2*iy)
            Ez = 0.5*(qez(2*ix,2*iy+1)+qez(2*ix+2,2*iy+1))
            Bx = qbx(2*ix+1,2*iy)
            By = qby(2*ix,2*iy+1)
            Bz = qbz(2*ix,2*iy)
            qMz(ix,iy) = (-1)*2*xi*c**2*(2*(Ex**2+Ey**2+Ez**2-c**2*(Bx**2+By**2+Bz**2))*Bz-7*(Ex*Bx+Ey*By+Ez*Bz)*Ez)
            ! write(*,*) Ex,Ey,Ez,Bx,By,Bz
         ENDDO

      ENDDO

      IF(x_min_boundary) THEN
         do iy=1-ng+1,ny+ng+1
            qMz(0,iy) = qMz(1,iy)
            qPy(1,iy) = 0.0_num
            qMy(0,iy) = qMy(1,iy)
            qPz(1,iy) = 0.0_num
         enddo
      ENDIF
      IF(x_max_boundary) THEN
         do iy=1-ng+1,ny+ng+1
            qMz(nx,iy) = qMz(nx-1,iy)
            qPy(nx,iy) = 0.0_num
            qMy(nx,iy) = qMy(nx-1,iy)
            qPz(nx,iy) = 0.0_num
         enddo
      ENDIF
   END SUBROUTINE calculate_PM_for_ex



   SUBROUTINE interpolation_for_ex0
      INTEGER :: ix,iy
      DO ix = 1-ng, nx+ng
         DO iy= 1-ng,ny+ng
            qex(2*ix,2*iy) = ex0(ix,iy)+static_ex(ix,iy)
            qey(2*ix,2*iy) = ey0(ix,iy)+static_ey(ix,iy)
            qez(2*ix,2*iy) = ez0(ix,iy)+static_ez(ix,iy)
            qbx(2*ix,2*iy) = bx0(ix,iy)+static_bx(ix,iy)
            qby(2*ix,2*iy) = by0(ix,iy)+static_by(ix,iy)
            qbz(2*ix,2*iy) = bz0(ix,iy)+static_bz(ix,iy)
         ENDDO
      ENDDO
      DO ix = 1-ng, nx+ng-1
         DO iy = 1-ng, ny+ng-1
            qex(2*ix+1,2*iy) = 0.5_num *(ex0(ix,iy)+ex0(ix+1,iy))+0.5_num * (static_ex(ix,iy)+ static_ex(ix+1,iy))
            qey(2*ix+1,2*iy) = 0.5_num *(ey0(ix,iy)+ey0(ix+1,iy))+0.5_num * (static_ey(ix,iy)+ static_ey(ix+1,iy))
            qez(2*ix+1,2*iy) = 0.5_num *(ez0(ix,iy)+ez0(ix+1,iy))+0.5_num * (static_ez(ix,iy)+ static_ez(ix+1,iy))
            qbx(2*ix+1,2*iy) = 0.5_num *(bx0(ix,iy)+bx0(ix+1,iy))+0.5_num * (static_bx(ix,iy)+ static_bx(ix+1,iy))
            qby(2*ix+1,2*iy) = 0.5_num *(by0(ix,iy)+by0(ix+1,iy))+0.5_num * (static_by(ix,iy)+ static_by(ix+1,iy))
            qbz(2*ix+1,2*iy) = 0.5_num *(bz0(ix,iy)+bz0(ix+1,iy))+0.5_num * (static_bz(ix,iy)+ static_bz(ix+1,iy))

            qex(2*ix,2*iy+1) = 0.5_num *(ex0(ix,iy)+ex0(ix,iy+1))+0.5_num * (static_ex(ix,iy)+ static_ex(ix,iy+1))
            qey(2*ix,2*iy+1) = 0.5_num *(ey0(ix,iy)+ey0(ix,iy+1))+0.5_num * (static_ey(ix,iy)+ static_ey(ix,iy+1))
            qez(2*ix,2*iy+1) = 0.5_num *(ez0(ix,iy)+ez0(ix,iy+1))+0.5_num * (static_ez(ix,iy)+ static_ez(ix,iy+1))
            qbx(2*ix,2*iy+1) = 0.5_num *(bx0(ix,iy)+bx0(ix,iy+1))+0.5_num * (static_bx(ix,iy)+ static_bx(ix,iy+1))
            qby(2*ix,2*iy+1) = 0.5_num *(by0(ix,iy)+by0(ix,iy+1))+0.5_num * (static_by(ix,iy)+ static_by(ix,iy+1))
            qbz(2*ix,2*iy+1) = 0.5_num *(bz0(ix,iy)+bz0(ix,iy+1))+0.5_num * (static_bz(ix,iy)+ static_bz(ix,iy+1))

            ! qex(2*ix+1,2*iy+1)=0.25_num*(ex0(ix,iy)+ex0(ix,iy+1)+ex0(ix+1,iy)+ex0(ix+1,iy+1))
            ! qey(2*ix+1,2*iy+1)=0.25_num*(ey0(ix,iy)+ey0(ix,iy+1)+ey0(ix+1,iy)+ey0(ix+1,iy+1))
            ! qez(2*ix+1,2*iy+1)=0.25_num*(ez0(ix,iy)+ez0(ix,iy+1)+ez0(ix+1,iy)+ez0(ix+1,iy+1))
            ! qbx(2*ix+1,2*iy+1)=0.25_num*(bx0(ix,iy)+bx0(ix,iy+1)+bx0(ix+1,iy)+bx0(ix+1,iy+1))
            ! qby(2*ix+1,2*iy+1)=0.25_num*(by0(ix,iy)+by0(ix,iy+1)+by0(ix+1,iy)+by0(ix+1,iy+1))
            ! qbz(2*ix+1,2*iy+1)=0.25_num*(bz0(ix,iy)+bz0(ix,iy+1)+bz0(ix+1,iy)+bz0(ix+1,iy+1))
         ENDDO
      ENDDO
      Do ix = 1-ng, nx+ng-1
         Do iy = 1-ng, ny+ng-1
             qex(2*ix+1,2*iy+1) =0.25_num * (qex(2*ix,2*iy+1)+qex(2*ix+2,2*iy+1)+qex(2*ix+1,2*iy)+qex(2*ix+1,2*iy+2))
             qey(2*ix+1,2*iy+1) =0.25_num * (qey(2*ix,2*iy+1)+qey(2*ix+2,2*iy+1)+qey(2*ix+1,2*iy)+qey(2*ix+1,2*iy+2))
             qez(2*ix+1,2*iy+1) =0.25_num * (qez(2*ix,2*iy+1)+qez(2*ix+2,2*iy+1)+qez(2*ix+1,2*iy)+qez(2*ix+1,2*iy+2))
             qbx(2*ix+1,2*iy+1) =0.25_num * (qbx(2*ix,2*iy+1)+qbx(2*ix+2,2*iy+1)+qbx(2*ix+1,2*iy)+qbx(2*ix+1,2*iy+2))
             qby(2*ix+1,2*iy+1) =0.25_num * (qby(2*ix,2*iy+1)+qby(2*ix+2,2*iy+1)+qby(2*ix+1,2*iy)+qby(2*ix+1,2*iy+2))
             qbz(2*ix+1,2*iy+1) =0.25_num * (qbz(2*ix,2*iy+1)+qbz(2*ix+2,2*iy+1)+qbz(2*ix+1,2*iy)+qbz(2*ix+1,2*iy+2))
         ENDDO
      ENDDO
   END SUBROUTINE interpolation_for_ex0

   SUBROUTINE calculate_PM_for_ex0
      INTEGER :: ix,iy
      REAL(num) :: xi
      REAL(num) :: Ex, Ey, Ez, Bx, By, Bz
      !xi = 2*alpha**2*epsilon0**2*h_bar**3/45/m0**4/c**5*1e8/2.285
      !xi = 5.71074693e-45_num
      !xi = 1.0e-15_num
      xi = ksi
      !write(*,*) xi
      DO ix = 1-ng+1, nx+ng-1
         DO iy = 1-ng+1, ny+ng-1
            !for Py, Mx
            !沿x方向做平均
            Ex = 0.5*(qex(2*ix,2*iy+1)+qex(2*ix-2,2*iy+1))
            Ey = qey(2*ix,2*iy)
            Ez = qez(2*ix,2*iy+1)
            Bx = qbx(2*ix,2*iy)
            By = 0.5*(qby(2*ix,2*iy+1)+qby(2*ix-2,2*iy+1))
            Bz = qbz(2*ix-1,2*iy)
            ! write(*,*) Ex,Ey,Ez,Bx,By,Bz
            qPy0(ix,iy) = 2*xi*(2*(Ex**2+Ey**2+Ez**2-c**2*(Bx**2+By**2+Bz**2))*Ey + 7*c**2*(Ex*Bx+Ey*By+Ez*Bz)*By)
            qMx0(ix,iy) = (-1)*2*xi*c**2*(2*(Ex**2+Ey**2+Ez**2-c**2*(Bx**2+By**2+Bz**2))*Bx-7*(Ex*Bx+Ey*By+Ez*Bz)*Ex)
            !write(*,*) 2*xi*(2*(Ex**2+Ey**2+Ez**2-c**2*(Bx**2+By**2+Bz**2))*Ey)
            !for Pz
            Ex = qex(2*ix-1,2*iy)
            Ey = qey(2*ix,2*iy-1)
            Ez = qez(2*ix,2*iy)
            Bx = qbx(2*ix,2*iy-1)
            By = qby(2*ix-1,2*iy)
            Bz = 0.5*(qbz(2*ix,2*iy-1)+qbz(2*ix-2,2*iy-1))
            ! write(*,*) Ex,Ey,Ez,Bx,By,Bz
            qPz0(ix,iy) = 2*xi*(2*(Ex**2+Ey**2+Ez**2-c**2*(Bx**2+By**2+Bz**2))*Ez + 7*c**2*(Ex*Bx+Ey*By+Ez*Bz)*Bz)
            !for Px, My
            Ex = qex(2*ix,2*iy)
            Ey = 0.5*(qey(2*ix,2*iy-1)+qey(2*ix+2,2*iy-1))
            Ez = qez(2*ix+1,2*iy)
            Bx = 0.5*(qbx(2*ix,2*iy-1)+qbx(2*ix+2,2*iy-1))
            By = qby(2*ix,2*iy)
            Bz = qbz(2*ix,2*iy-1)
            qPx0(ix,iy) = 2*xi*(2*(Ex**2+Ey**2+Ez**2-c**2*(Bx**2+By**2+Bz**2))*Ex + 7*c**2*(Ex*Bx+Ey*By+Ez*Bz)*Bx)
            qMy0(ix,iy) = (-1)*2*xi*c**2*(2*(Ex**2+Ey**2+Ez**2-c**2*(Bx**2+By**2+Bz**2))*By-7*(Ex*Bx+Ey*By+Ez*Bz)*Ey)
            ! write(*,*) Ex,Ey,Ez,Bx,By,Bz
            !write(*,*) qPx0(ix),qPy0(ix),qPz0(ix)
            !write(*,*) qMx0(ix),qMy0(ix),qMz0(ix)
            !for Mz
            Ex = qex(2*ix,2*iy+1)
            Ey = qey(2*ix+1,2*iy)
            Ez = 0.5*(qez(2*ix,2*iy+1)+qez(2*ix+2,2*iy+1))
            Bx = qbx(2*ix+1,2*iy)
            By = qby(2*ix,2*iy+1)
            Bz = qbz(2*ix,2*iy)
            qMz0(ix,iy) = (-1)*2*xi*c**2*(2*(Ex**2+Ey**2+Ez**2-c**2*(Bx**2+By**2+Bz**2))*Bz-7*(Ex*Bx+Ey*By+Ez*Bz)*Ez)
            ! write(*,*) Ex,Ey,Ez,Bx,By,Bz
         ENDDO

      ENDDO
      IF(x_min_boundary) THEN
         do iy=1-ng+1,ny+ng+1
            qMz0(0,iy) = qMz0(1,iy)
            qPy0(1,iy) = 0.0_num
            qMy0(0,iy) = qMy0(1,iy)
            qPz0(1,iy) = 0.0_num
         enddo
      ENDIF
      IF(x_max_boundary) THEN
         do iy=1-ng+1,ny+ng+1
            qMz0(nx,iy) = qMz0(nx-1,iy)
            qPy0(nx,iy) = 0.0_num
            qMy0(nx,iy) = qMy0(nx-1,iy)
            qPz0(nx,iy) = 0.0_num
         enddo
      ENDIF
   END SUBROUTINE calculate_PM_for_ex0

   SUBROUTINE update_e_field_zero
      INTEGER :: ix,iy
      REAL(num) :: cx1
      REAL(num) :: cy1
      cx1 = cnx
      cy1 = cny
      !场二阶推进
      DO iy = 0, ny
         DO ix = 0, nx
            ex0(ix, iy) = ex(ix, iy) &
               + cy1 * (bz(ix  , iy  ) - bz(ix  , iy-1)) &
               - fac * jx(ix, iy)

            ey0(ix, iy) = ey(ix, iy) &
               - cx1 * (bz(ix  , iy  ) - bz(ix-1, iy  )) &
               - fac * jy(ix, iy)

            ez0(ix, iy) = ez(ix, iy) &
               + cx1 * (by(ix  , iy  ) - by(ix-1, iy  )) &
               - cy1 * (bx(ix  , iy  ) - bx(ix  , iy-1)) &
               - fac * jz(ix, iy)
         END DO
      END DO
   END SUBROUTINE update_e_field_zero

   SUBROUTINE update_b_field_zero
      !=========
      INTEGER :: ix,iy
      REAL(num) :: cx1
      REAL(num) :: cy1
      cx1 = hdtx
      cy1 = hdty
      DO iy = 0, ny
         DO ix = 0, nx
            bx0(ix, iy) = bx(ix, iy) &
               - cy1 * (ez0(ix  , iy+1) - ez0(ix  , iy  ))

            by0(ix, iy) = by(ix, iy) &
               + cx1 * (ez0(ix+1, iy  ) - ez0(ix  , iy  ))

            bz0(ix, iy) = bz(ix, iy) &
               - cx1 * (ey0(ix+1, iy  ) - ey0(ix  , iy  )) &
               + cy1 * (ex0(ix  , iy+1) - ex0(ix  , iy  ))
         END DO
      END DO

   END SUBROUTINE update_b_field_zero

   SUBROUTINE update_e_field_vcm
      !场二阶推进
      INTEGER :: ix,iy
      REAL(num) :: cx1
      REAL(num) :: cy1
      cx1 = cnx
      cy1 = cny
      !场二阶推进
      DO iy = 0, ny
         DO ix = 0, nx
            ex1(ix, iy) = ex(ix, iy) &
               + cy1 * ( (bz(ix  , iy  ) - bz(ix  , iy-1)) &
               -mu0*(qMz0(ix,iy)-qMz0(ix,iy-1) ) ) &
               - fac * jx(ix, iy)&
               - ( qPx0(ix,iy) - qPx(ix,iy) )/epsilon0

            ey1(ix, iy) = ey(ix, iy) &
               - cx1 *( (bz(ix  , iy  ) - bz(ix-1, iy  )) &
               -mu0*(qMz0(ix,iy)-qMz0(ix-1,iy)) ) &
               - fac * jy(ix, iy)&
               - ( qPy0(ix,iy) - qPy(ix,iy) )/epsilon0
            ez1(ix, iy) = ez(ix, iy) &
               + cx1 *( (by(ix  , iy  ) - by(ix-1, iy  )) &
               -mu0*(qMy0(ix  , iy) - qMy0(ix-1  , iy) ) )&
               - cy1 * ( (bx(ix  , iy  ) - bx(ix  , iy-1)) &
               -mu0*(qMx0(ix  , iy) -qMx0(ix  , iy-1) ) )&
               - fac * jz(ix, iy)&
               - ( qPz0(ix,iy) - qPz(ix,iy) )/epsilon0
            ! write(*,*)ix
            ! write(*,*)iy
            !delta_ex(ix,iy)=ex1(ix, iy)-ex(ix,iy)
            !delta_ey(ix,iy)=ey1(ix, iy)-ey(ix,iy)
            !delta_ez(ix,iy)=ez1(ix, iy)-ez(ix,iy)
         END DO
      END DO
   END SUBROUTINE update_e_field_vcm


   SUBROUTINE update_b_field_vcm
      INTEGER :: ix,iy
      REAL(num) :: cx1
      REAL(num) :: cy1
      cx1 = hdtx
      cy1 = hdty

      DO iy = 0, ny
         DO ix = 0, nx
            bx1(ix, iy) = bx(ix, iy) &
               - cy1 * (ez1(ix  , iy+1) - ez1(ix  , iy  ))

            by1(ix, iy) = by(ix, iy) &
               + cx1 * (ez1(ix+1, iy  ) - ez1(ix  , iy  ))

            bz1(ix, iy) = bz(ix, iy) &
               - cx1 * (ey1(ix+1, iy  ) - ey1(ix  , iy  )) &
               + cy1 * (ex1(ix  , iy+1) - ex1(ix  , iy  ))
         END DO
      END DO

   END SUBROUTINE update_b_field_vcm


!=============================================
!============非vcm解法========================
!=============================================
   SUBROUTINE set_field_order(order)

      INTEGER, INTENT(IN) :: order

      field_order = order
      fng = field_order / 2

      IF (field_order == 2) THEN
         cfl = 1.0_num
      ELSE IF (field_order == 4) THEN
         cfl = 6.0_num / 7.0_num
      ELSE
         cfl = 120.0_num / 149.0_num
      END IF

   END SUBROUTINE set_field_order



   SUBROUTINE set_maxwell_solver

      REAL(num) :: delta, dx_cdt

      IF (maxwell_solver == c_maxwell_solver_custom) THEN
         alphax = 1.0_num - 2.0_num * betaxy - 3.0_num * deltax
         alphay = 1.0_num - 2.0_num * betayx - 3.0_num * deltay

      ELSE IF (maxwell_solver == c_maxwell_solver_lehe_x) THEN
         ! R. Lehe et al., Phys. Rev. ST Accel. Beams 16, 021301 (2013)
         dx_cdt = dx / (c * dt)
         betaxy = 0.125_num * (dx / dy)**2
         betayx = 0.125_num
         deltax = 0.25_num * (1.0_num - dx_cdt**2 * SIN(0.5_num * pi / dx_cdt)**2)
         deltay = 0.0_num
         alphax = 1.0_num - 2.0_num * betaxy - 3.0_num * deltax
         alphay = 1.0_num - 2.0_num * betayx

      ELSE IF (maxwell_solver == c_maxwell_solver_lehe_y) THEN
         dx_cdt = dy / (c * dt)
         betayx = 0.125_num * (dy / dx)**2
         betaxy = 0.125_num
         deltax = 0.0_num
         deltay = 0.25_num * (1.0_num - dx_cdt**2 * SIN(0.5_num * pi / dx_cdt)**2)
         alphax = 1.0_num - 2.0_num * betaxy
         alphay = 1.0_num - 2.0_num * betayx - 3.0_num * deltay

      ELSE IF (maxwell_solver == c_maxwell_solver_pukhov) THEN
         ! A. Pukhov, Journal of Plasma Physics 61, 425-433 (1999)
         delta = MIN(dx, dy)

         betayx = 0.125_num * (delta / dx)**2
         betaxy = 0.125_num * (delta / dy)**2
         deltax = 0.0_num
         deltay = 0.0_num
         alphax = 1.0_num - 2.0_num * betaxy
         alphay = 1.0_num - 2.0_num * betayx
      END IF

      IF (rank == 0 .AND. maxwell_solver /= c_maxwell_solver_yee) THEN
         PRINT*, 'Maxwell solver set to the following parameters:'
         PRINT'(A9, 2F14.9)', 'alpha =', alphax, alphay
         PRINT'(A9, 1F14.9)', 'betax =', betaxy
         PRINT'(A9, 1F14.9)', 'betay =', betayx
         PRINT'(A9, 2F14.9)', 'delta =', deltax, deltay
         PRINT'(A9, 1F14.9)', 'c*dt/dx = ', dt * c / dx
         PRINT*
      END IF

   END SUBROUTINE set_maxwell_solver



   SUBROUTINE update_e_field

      INTEGER :: ix, iy
      REAL(num) :: cpml_x, cpml_y
      REAL(num) :: c1, c2, c3
      REAL(num) :: cx1, cx2, cx3
      REAL(num) :: cy1, cy2, cy3

      IF (cpml_boundaries) THEN
         IF (field_order == 2) THEN
            DO iy = 0, ny
               cy1 = cny / cpml_kappa_ey(iy)
               DO ix = 0, nx
                  cx1 = cnx / cpml_kappa_ex(ix)

                  ex(ix, iy) = ex(ix, iy) &
                     + cy1 * (bz(ix  , iy  ) - bz(ix  , iy-1)) &
                     - fac * jx(ix, iy)

                  ey(ix, iy) = ey(ix, iy) &
                     - cx1 * (bz(ix  , iy  ) - bz(ix-1, iy  )) &
                     - fac * jy(ix, iy)

                  ez(ix, iy) = ez(ix, iy) &
                     + cx1 * (by(ix  , iy  ) - by(ix-1, iy  )) &
                     - cy1 * (bx(ix  , iy  ) - bx(ix  , iy-1)) &
                     - fac * jz(ix, iy)
               END DO
            END DO
         ELSE IF (field_order == 4) THEN
            c1 = 9.0_num / 8.0_num
            c2 = -1.0_num / 24.0_num

            DO iy = 0, ny
               cpml_y = cny / cpml_kappa_ey(iy)
               cy1 = c1 * cpml_y
               cy2 = c2 * cpml_y
               DO ix = 0, nx
                  cpml_x = cnx / cpml_kappa_ex(ix)
                  cx1 = c1 * cpml_x
                  cx2 = c2 * cpml_x

                  ex(ix, iy) = ex(ix, iy) &
                     + cy1 * (bz(ix  , iy  ) - bz(ix  , iy-1)) &
                     + cy2 * (bz(ix  , iy+1) - bz(ix  , iy-2)) &
                     - fac * jx(ix, iy)

                  ey(ix, iy) = ey(ix, iy) &
                     - cx1 * (bz(ix  , iy  ) - bz(ix-1, iy  )) &
                     - cx2 * (bz(ix+1, iy  ) - bz(ix-2, iy  )) &
                     - fac * jy(ix, iy)

                  ez(ix, iy) = ez(ix, iy) &
                     + cx1 * (by(ix  , iy  ) - by(ix-1, iy  )) &
                     + cx2 * (by(ix+1, iy  ) - by(ix-2, iy  )) &
                     - cy1 * (bx(ix  , iy  ) - bx(ix  , iy-1)) &
                     - cy2 * (bx(ix  , iy+1) - bx(ix  , iy-2)) &
                     - fac * jz(ix, iy)
               END DO
            END DO
         ELSE
            c1 = 75.0_num / 64.0_num
            c2 = -25.0_num / 384.0_num
            c3 = 3.0_num / 640.0_num

            DO iy = 0, ny
               cpml_y = cny / cpml_kappa_ey(iy)
               cy1 = c1 * cpml_y
               cy2 = c2 * cpml_y
               cy3 = c3 * cpml_y
               DO ix = 0, nx
                  cpml_x = cnx / cpml_kappa_ex(ix)
                  cx1 = c1 * cpml_x
                  cx2 = c2 * cpml_x
                  cx3 = c3 * cpml_x

                  ex(ix, iy) = ex(ix, iy) &
                     + cy1 * (bz(ix  , iy  ) - bz(ix  , iy-1)) &
                     + cy2 * (bz(ix  , iy+1) - bz(ix  , iy-2)) &
                     + cy3 * (bz(ix  , iy+2) - bz(ix  , iy-3)) &
                     - fac * jx(ix, iy)

                  ey(ix, iy) = ey(ix, iy) &
                     - cx1 * (bz(ix  , iy  ) - bz(ix-1, iy  )) &
                     - cx2 * (bz(ix+1, iy  ) - bz(ix-2, iy  )) &
                     - cx3 * (bz(ix+2, iy  ) - bz(ix-3, iy  )) &
                     - fac * jy(ix, iy)

                  ez(ix, iy) = ez(ix, iy) &
                     + cx1 * (by(ix  , iy  ) - by(ix-1, iy  )) &
                     + cx2 * (by(ix+1, iy  ) - by(ix-2, iy  )) &
                     + cx3 * (by(ix+2, iy  ) - by(ix-3, iy  )) &
                     - cy1 * (bx(ix  , iy  ) - bx(ix  , iy-1)) &
                     - cy2 * (bx(ix  , iy+1) - bx(ix  , iy-2)) &
                     - cy3 * (bx(ix  , iy+2) - bx(ix  , iy-3)) &
                     - fac * jz(ix, iy)
               END DO
            END DO
         END IF

         CALL cpml_advance_e_currents(hdt)
      ELSE
         IF (field_order == 2) THEN
            cx1 = cnx
            cy1 = cny

            DO iy = 0, ny
               DO ix = 0, nx
                  ex(ix, iy) = ex(ix, iy) &
                     + cy1 * (bz(ix  , iy  ) - bz(ix  , iy-1)) &
                     - fac * jx(ix, iy)

                  ey(ix, iy) = ey(ix, iy) &
                     - cx1 * (bz(ix  , iy  ) - bz(ix-1, iy  )) &
                     - fac * jy(ix, iy)

                  ez(ix, iy) = ez(ix, iy) &
                     + cx1 * (by(ix  , iy  ) - by(ix-1, iy  )) &
                     - cy1 * (bx(ix  , iy  ) - bx(ix  , iy-1)) &
                     - fac * jz(ix, iy)
               END DO
            END DO
         ELSE IF (field_order == 4) THEN
            c1 = 9.0_num / 8.0_num
            c2 = -1.0_num / 24.0_num

            cx1 = c1 * cnx
            cx2 = c2 * cnx
            cy1 = c1 * cny
            cy2 = c2 * cny

            DO iy = 0, ny
               DO ix = 0, nx
                  ex(ix, iy) = ex(ix, iy) &
                     + cy1 * (bz(ix  , iy  ) - bz(ix  , iy-1)) &
                     + cy2 * (bz(ix  , iy+1) - bz(ix  , iy-2)) &
                     - fac * jx(ix, iy)

                  ey(ix, iy) = ey(ix, iy) &
                     - cx1 * (bz(ix  , iy  ) - bz(ix-1, iy  )) &
                     - cx2 * (bz(ix+1, iy  ) - bz(ix-2, iy  )) &
                     - fac * jy(ix, iy)

                  ez(ix, iy) = ez(ix, iy) &
                     + cx1 * (by(ix  , iy  ) - by(ix-1, iy  )) &
                     + cx2 * (by(ix+1, iy  ) - by(ix-2, iy  )) &
                     - cy1 * (bx(ix  , iy  ) - bx(ix  , iy-1)) &
                     - cy2 * (bx(ix  , iy+1) - bx(ix  , iy-2)) &
                     - fac * jz(ix, iy)
               END DO
            END DO
         ELSE
            c1 = 75.0_num / 64.0_num
            c2 = -25.0_num / 384.0_num
            c3 = 3.0_num / 640.0_num

            cx1 = c1 * cnx
            cx2 = c2 * cnx
            cx3 = c3 * cnx
            cy1 = c1 * cny
            cy2 = c2 * cny
            cy3 = c3 * cny

            DO iy = 0, ny
               DO ix = 0, nx
                  ex(ix, iy) = ex(ix, iy) &
                     + cy1 * (bz(ix  , iy  ) - bz(ix  , iy-1)) &
                     + cy2 * (bz(ix  , iy+1) - bz(ix  , iy-2)) &
                     + cy3 * (bz(ix  , iy+2) - bz(ix  , iy-3)) &
                     - fac * jx(ix, iy)

                  ey(ix, iy) = ey(ix, iy) &
                     - cx1 * (bz(ix  , iy  ) - bz(ix-1, iy  )) &
                     - cx2 * (bz(ix+1, iy  ) - bz(ix-2, iy  )) &
                     - cx3 * (bz(ix+2, iy  ) - bz(ix-3, iy  )) &
                     - fac * jy(ix, iy)

                  ez(ix, iy) = ez(ix, iy) &
                     + cx1 * (by(ix  , iy  ) - by(ix-1, iy  )) &
                     + cx2 * (by(ix+1, iy  ) - by(ix-2, iy  )) &
                     + cx3 * (by(ix+2, iy  ) - by(ix-3, iy  )) &
                     - cy1 * (bx(ix  , iy  ) - bx(ix  , iy-1)) &
                     - cy2 * (bx(ix  , iy+1) - bx(ix  , iy-2)) &
                     - cy3 * (bx(ix  , iy+2) - bx(ix  , iy-3)) &
                     - fac * jz(ix, iy)
               END DO
            END DO
         END IF
      END IF

   END SUBROUTINE update_e_field



   SUBROUTINE update_b_field

      INTEGER :: ix, iy
      REAL(num) :: cpml_x, cpml_y
      REAL(num) :: c1, c2, c3
      REAL(num) :: cx1, cx2, cx3
      REAL(num) :: cy1, cy2, cy3

      IF (cpml_boundaries) THEN
         IF (field_order == 2) THEN
            IF (maxwell_solver == c_maxwell_solver_yee) THEN
               DO iy = 0, ny
                  cy1 = hdty / cpml_kappa_by(iy)
                  DO ix = 0, nx
                     cx1 = hdtx / cpml_kappa_bx(ix)

                     bx(ix, iy) = bx(ix, iy) &
                        - cy1 * (ez(ix  , iy+1) - ez(ix  , iy  ))

                     by(ix, iy) = by(ix, iy) &
                        + cx1 * (ez(ix+1, iy  ) - ez(ix  , iy  ))

                     bz(ix, iy) = bz(ix, iy) &
                        - cx1 * (ey(ix+1, iy  ) - ey(ix  , iy  )) &
                        + cy1 * (ex(ix  , iy+1) - ex(ix  , iy  ))
                  END DO
               END DO
            ELSE
               DO iy = 0, ny
                  cy1 = hdty / cpml_kappa_by(iy)
                  DO ix = 0, nx
                     cx1 = hdtx / cpml_kappa_bx(ix)

                     bx(ix, iy) = bx(ix, iy) &
                        - cy1 * (alphay * (ez(ix  , iy+1) - ez(ix  , iy  ))  &
                        + betayx * (ez(ix+1, iy+1) - ez(ix+1, iy  )   &
                        + ez(ix-1, iy+1) - ez(ix-1, iy  ))  &
                        + deltay * (ez(ix  , iy+2) - ez(ix  , iy-1)))

                     by(ix, iy) = by(ix, iy) &
                        + cx1 * (alphax * (ez(ix+1, iy  ) - ez(ix  , iy  ))  &
                        + betaxy * (ez(ix+1, iy+1) - ez(ix  , iy+1)   &
                        + ez(ix+1, iy-1) - ez(ix  , iy-1))  &
                        + deltax * (ez(ix+2, iy  ) - ez(ix-1, iy  )))

                     bz(ix, iy) = bz(ix, iy) &
                        - cx1 * (alphax * (ey(ix+1, iy  ) - ey(ix  , iy  ))  &
                        + betaxy * (ey(ix+1, iy+1) - ey(ix  , iy+1)   &
                        + ey(ix+1, iy-1) - ey(ix  , iy-1))  &
                        + deltax * (ey(ix+2, iy  ) - ey(ix-1, iy  ))) &
                        + cy1 * (alphay * (ex(ix  , iy+1) - ex(ix  , iy  ))  &
                        + betayx * (ex(ix+1, iy+1) - ex(ix+1, iy  )   &
                        + ex(ix-1, iy+1) - ex(ix-1, iy  ))  &
                        + deltay * (ex(ix  , iy+2) - ex(ix  , iy-1)))
                  END DO
               END DO
            END IF
         ELSE IF (field_order == 4) THEN
            c1 = 9.0_num / 8.0_num
            c2 = -1.0_num / 24.0_num

            DO iy = 0, ny
               cpml_y = hdty / cpml_kappa_by(iy)
               cy1 = c1 * cpml_y
               cy2 = c2 * cpml_y
               DO ix = 0, nx
                  cpml_x = hdtx / cpml_kappa_bx(ix)
                  cx1 = c1 * cpml_x
                  cx2 = c2 * cpml_x

                  bx(ix, iy) = bx(ix, iy) &
                     - cy1 * (ez(ix  , iy+1) - ez(ix  , iy  )) &
                     - cy2 * (ez(ix  , iy+2) - ez(ix  , iy-1))

                  by(ix, iy) = by(ix, iy) &
                     + cx1 * (ez(ix+1, iy  ) - ez(ix  , iy  )) &
                     + cx2 * (ez(ix+2, iy  ) - ez(ix-1, iy  ))

                  bz(ix, iy) = bz(ix, iy) &
                     - cx1 * (ey(ix+1, iy  ) - ey(ix  , iy  )) &
                     - cx2 * (ey(ix+2, iy  ) - ey(ix-1, iy  )) &
                     + cy1 * (ex(ix  , iy+1) - ex(ix  , iy  )) &
                     + cy2 * (ex(ix  , iy+2) - ex(ix  , iy-1))
               END DO
            END DO
         ELSE
            c1 = 75.0_num / 64.0_num
            c2 = -25.0_num / 384.0_num
            c3 = 3.0_num / 640.0_num

            DO iy = 0, ny
               cpml_y = hdty / cpml_kappa_by(iy)
               cy1 = c1 * cpml_y
               cy2 = c2 * cpml_y
               cy3 = c3 * cpml_y
               DO ix = 0, nx
                  cpml_x = hdtx / cpml_kappa_bx(ix)
                  cx1 = c1 * cpml_x
                  cx2 = c2 * cpml_x
                  cx3 = c3 * cpml_x

                  bx(ix, iy) = bx(ix, iy) &
                     - cy1 * (ez(ix  , iy+1) - ez(ix  , iy  )) &
                     - cy2 * (ez(ix  , iy+2) - ez(ix  , iy-1)) &
                     - cy3 * (ez(ix  , iy+3) - ez(ix  , iy-2))

                  by(ix, iy) = by(ix, iy) &
                     + cx1 * (ez(ix+1, iy  ) - ez(ix  , iy  )) &
                     + cx2 * (ez(ix+2, iy  ) - ez(ix-1, iy  )) &
                     + cx3 * (ez(ix+3, iy  ) - ez(ix-2, iy  ))

                  bz(ix, iy) = bz(ix, iy) &
                     - cx1 * (ey(ix+1, iy  ) - ey(ix  , iy  )) &
                     - cx2 * (ey(ix+2, iy  ) - ey(ix-1, iy  )) &
                     - cx3 * (ey(ix+3, iy  ) - ey(ix-2, iy  )) &
                     + cy1 * (ex(ix  , iy+1) - ex(ix  , iy  )) &
                     + cy2 * (ex(ix  , iy+2) - ex(ix  , iy-1)) &
                     + cy3 * (ex(ix  , iy+3) - ex(ix  , iy-2))
               END DO
            END DO
         END IF

         CALL cpml_advance_b_currents(hdt)
      ELSE
         IF (field_order == 2) THEN
            cx1 = hdtx
            cy1 = hdty

            IF (maxwell_solver == c_maxwell_solver_yee) THEN
               DO iy = 0, ny
                  DO ix = 0, nx
                     bx(ix, iy) = bx(ix, iy) &
                        - cy1 * (ez(ix  , iy+1) - ez(ix  , iy  ))

                     by(ix, iy) = by(ix, iy) &
                        + cx1 * (ez(ix+1, iy  ) - ez(ix  , iy  ))

                     bz(ix, iy) = bz(ix, iy) &
                        - cx1 * (ey(ix+1, iy  ) - ey(ix  , iy  )) &
                        + cy1 * (ex(ix  , iy+1) - ex(ix  , iy  ))
                  END DO
               END DO
            ELSE
               DO iy = 0, ny
                  DO ix = 0, nx
                     bx(ix, iy) = bx(ix, iy) &
                        - cy1 * (alphay * (ez(ix  , iy+1) - ez(ix  , iy  ))  &
                        + betayx * (ez(ix+1, iy+1) - ez(ix+1, iy  )   &
                        + ez(ix-1, iy+1) - ez(ix-1, iy  ))  &
                        + deltay * (ez(ix  , iy+2) - ez(ix  , iy-1)))

                     by(ix, iy) = by(ix, iy) &
                        + cx1 * (alphax * (ez(ix+1, iy  ) - ez(ix  , iy  ))  &
                        + betaxy * (ez(ix+1, iy+1) - ez(ix  , iy+1)   &
                        + ez(ix+1, iy-1) - ez(ix  , iy-1))  &
                        + deltax * (ez(ix+2, iy  ) - ez(ix-1, iy  )))

                     bz(ix, iy) = bz(ix, iy) &
                        - cx1 * (alphax * (ey(ix+1, iy  ) - ey(ix  , iy  ))  &
                        + betaxy * (ey(ix+1, iy+1) - ey(ix  , iy+1)   &
                        + ey(ix+1, iy-1) - ey(ix  , iy-1))  &
                        + deltax * (ey(ix+2, iy  ) - ey(ix-1, iy  ))) &
                        + cy1 * (alphay * (ex(ix  , iy+1) - ex(ix  , iy  ))  &
                        + betayx * (ex(ix+1, iy+1) - ex(ix+1, iy  )   &
                        + ex(ix-1, iy+1) - ex(ix-1, iy  ))  &
                        + deltay * (ex(ix  , iy+2) - ex(ix  , iy-1)))
                  END DO
               END DO
            END IF
         ELSE IF (field_order == 4) THEN
            c1 = 9.0_num / 8.0_num
            c2 = -1.0_num / 24.0_num

            cx1 = c1 * hdtx
            cx2 = c2 * hdtx
            cy1 = c1 * hdty
            cy2 = c2 * hdty

            DO iy = 0, ny
               DO ix = 0, nx
                  bx(ix, iy) = bx(ix, iy) &
                     - cy1 * (ez(ix  , iy+1) - ez(ix  , iy  )) &
                     - cy2 * (ez(ix  , iy+2) - ez(ix  , iy-1))

                  by(ix, iy) = by(ix, iy) &
                     + cx1 * (ez(ix+1, iy  ) - ez(ix  , iy  )) &
                     + cx2 * (ez(ix+2, iy  ) - ez(ix-1, iy  ))

                  bz(ix, iy) = bz(ix, iy) &
                     - cx1 * (ey(ix+1, iy  ) - ey(ix  , iy  )) &
                     - cx2 * (ey(ix+2, iy  ) - ey(ix-1, iy  )) &
                     + cy1 * (ex(ix  , iy+1) - ex(ix  , iy  )) &
                     + cy2 * (ex(ix  , iy+2) - ex(ix  , iy-1))
               END DO
            END DO
         ELSE
            c1 = 75.0_num / 64.0_num
            c2 = -25.0_num / 384.0_num
            c3 = 3.0_num / 640.0_num

            cx1 = c1 * hdtx
            cx2 = c2 * hdtx
            cx3 = c3 * hdtx
            cy1 = c1 * hdty
            cy2 = c2 * hdty
            cy3 = c3 * hdty

            DO iy = 0, ny
               DO ix = 0, nx
                  bx(ix, iy) = bx(ix, iy) &
                     - cy1 * (ez(ix  , iy+1) - ez(ix  , iy  )) &
                     - cy2 * (ez(ix  , iy+2) - ez(ix  , iy-1)) &
                     - cy3 * (ez(ix  , iy+3) - ez(ix  , iy-2))

                  by(ix, iy) = by(ix, iy) &
                     + cx1 * (ez(ix+1, iy  ) - ez(ix  , iy  )) &
                     + cx2 * (ez(ix+2, iy  ) - ez(ix-1, iy  )) &
                     + cx3 * (ez(ix+3, iy  ) - ez(ix-2, iy  ))

                  bz(ix, iy) = bz(ix, iy) &
                     - cx1 * (ey(ix+1, iy  ) - ey(ix  , iy  )) &
                     - cx2 * (ey(ix+2, iy  ) - ey(ix-1, iy  )) &
                     - cx3 * (ey(ix+3, iy  ) - ey(ix-2, iy  )) &
                     + cy1 * (ex(ix  , iy+1) - ex(ix  , iy  )) &
                     + cy2 * (ex(ix  , iy+2) - ex(ix  , iy-1)) &
                     + cy3 * (ex(ix  , iy+3) - ex(ix  , iy-2))
               END DO
            END DO
         END IF
      END IF

   END SUBROUTINE update_b_field


END MODULE fields
